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The dynamics of particles in turbulence when the particle-size is larger than the dissipa- 
tive scale of the carrier flow is studied. Recent experiments have highlighted signatures 
of particles finiteness on their statistical properties, namely a decrease of their acceler- 
ation variance, an increase of correlation times -at increasing the particles size- and an 
independence of the probability density function of the acceleration once normalized to 
their variance. These effects are not captured by point particle models. By means of a 
detailed comparison between numerical simulations and experimental data, we show that 
a more accurate model is obtained once Faxen corrections are included. 



The study of Lagrangian turbulence and of t urbulent transport of material particles 
has received growing interest in recent years (Toschi & Bodcnschatz (2009)). Modern 
experimental techniques (based on synchronization of multiple fast cameras, or ultra- 
sonic/laser Do££lB£vel£cmietr^ allow nowadays t o fully resolve particle tracks in turbu- 
lent flows (iLa Porta et oil (l200l|) ; [Mordant et al} (|200l[ ): Xu et al. (2006); Berg (2006h: 



Volk et al\ (|2008a )). These techniques have opened the way toward a systematic study 



of the dynamics of material (or inertial) particles. When the particle density is dif- 
ferent from the one of the carri er fluid, a rich p t ienorn e nology emerges, s u ch as par- 
ticle clustering and segregation ( Squires fc EatonI (jl99lh : ICalzavarini et al\ (|2008&l |ca)). 
Numerical studies have proven to be essential tools in complementing and benchmark- 
ing early days experimental data: investi gations of fluid tracer s dynamics have show n 



remarkable agreeme nt with experiments (|Mordant et al\ (|2004[ ) ; I Arneodo et al\ (|2008l ): 



tg reeme ni 

Biferale et al\ (l2008l) ). Lagrangian numerical studies through Direct Numerical Simula- 



tions (DNS) of very small - computationally assumed to be pointwise - particles have 
also shown encouragi n g consistency with exper i menta l measurements for in e rtial parti- 
cles ('Bee et o/J ( 2006 ): Avvalasomavaiula et al. (2006); Salazar et al. ( 2008 ): Volk et al 



(2008 a)). However, in many situations the size of the particles is not small with respect 
to turbulence scales, in particular the dissipative scale rj. One example is Plankton which, 
while neutrally buoyant, cannot be considered as a tracer because of its size in the order 
of few dissipative scales. Typical marine and atmospherical environmental flows have 
ri ~ O(10)/xm. 

The statistics of particle accelerations, which directly reflects the action of hydrody- 
namical forces, has been used to experimentally assess the statistical signature of "large" 
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spherical particles, i.e. whose d iamet er D is larger t h an th e smallest turbulence scale 
■q. Recent studies ( Voth et al. (20o2); Qureshi et al. (l2007h ) and d etailed comparison 
between experiments and numerical simulations ( Volk etaU (2008a)) have shown that 
finite-sized neutrally-buoyant particles cannot be modeled as pointwise in numerical stud- 
ies. Features which have been clearly associated with a finite particle size are: 

i) For neutrally buoyant particles with D > rj the acceleration variance decreases at 
increasing the particle size. A scaling law behavior, (a^) ~ e^/'^Z)"^/'^ (with e being the 
energy dissipation rate), has been suggeste d on the basis of Kolmogorov 1941 (K41) tur- 
bulence phenomenology ( Voth et al. ( 20021 ): [Oureshi et al\ ( 2007l )l. 

ii) The normalized acceleration probability density function (PDF) depends at best very 
weakly o n the particl e diam e ter. Its shape can be fitted with stretched exponential func- 
tions (see lVoth et~d\ (|2002l ): lOureshi et al\ (|2007l )). 

iii) The autocorrelatio n function o f acceleration shows increased correlation time with 
increasing particle-size (jVolk et all ( ,2008 a)). 

While experimentally it is easier to study large {D > rj) particles, theoretically (and 
therefore computationally) this turns out a far more difficult task. Our aim in this article 
is to study the novel features associated with finite particles size in developed turbu- 
lent flows, while presenting an improved numerical model cap able to solve most of the 
discrepancies between experiments and simulations noticed in ( Volk et~al\ ( 2008qI )1. We 
show that qualitatively and quantitatively the new features are well captured by an equa- 
tion of motion which takes into account the effect of the nonuniformity of the flow at 
the particle-scale. To our knowledge the impact on acceleration statistics of such forces, 
known since a long time as Faxen corrections ({Faxen (1922)) has never been considered. 

The article is organized as follows. First we comment on the problems of obtaining 
an equation of motion for finite-sized particles. We examine the approximation on which 
point-particles equations rely and discuss two highly simplified models for the dynam- 
ics of small {D < rj) and finite-sized {D > rj) particles. Section [2] gives the numerical 
implementation of the proposed Faxen-corrected model. In section |3] we show basic phys- 
ical differences between the statistics of particle acceleration given by numerics with or 
without Faxen corrections. Section |3] contains the comparison of the model against exper- 
imental results, focusing on neutrally buoyant particles. Finally in Sec. [5] we summarize 
the results, we critically review the model and discuss how it can be improved. 



1. Equation of motion for finite-sized particle in turbulence 

Many studies on fine particulate flows have based particle's description on an equation 
- referred to as Maxey-Riley-Gatignol - which is an exact derivation of the forces on 
a particle in a nonuniform unsteady flow in the limit of vanishing Reynolds numbers 
Rep = Dvs/v and Res — D^T/i', where Vs is the sli p particle veloci t y resp e ct to the 



fluid a nd F = |Vu| the typical shear-scale in the flow (jMaxev fc RilevI |l983h : iGatignol 
( 19831 )'). In the net hydrodynamical force acting on a particle given by this equation 
one recognizes several contributions: the steady Stokes drag, the fluid acceleration force 
(sum of the pressure gradient and the dissipative forces on the fluid), the added mass, the 
buoyancy, the history Basset-Boussinesq force, and Faxen corrections. When the control 
parameters Rcp and Res become finite, the non-linearity of the f low dynamics in t he 
vicinity of the particle must be taken into account (see the review bv lMichaelide^ HM^)). 
A n expression for the a dded mass term which is correct at any Rep value has been derived 



bv iAuton et al. ( 1988[ ). But much more complicated is the situation for the other forces 



involved. The drag term becomes Rep depende nt and emp irical expressions based on 
numerical computations have been proposed (see Clift et al. (|l978) ). Furthermore, a lift 
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force appears at finite values of Rcp and Res- This force is notably hard to model because 
of the non-linear combination of shear and vorticity, and approximate expressions based 
on Saffman (small Rcp) and Lighthill-Auton (l arge Re^)) mechanisms are often used in 
studies (see e.g. discussion on lift on bubbles by iMagnaudet fc Legendr 



Theoretical and numerical studies of fine disperse multi-phase flows, which aim at de- 
scribing the behavior of a large number of particles, have adopted simplifi ed models where 
the su b -dom i nant t erms in Maxey-Riley-Gatignol equation are neglected te alkovsky et al 



( 200l[ ): iBed ( 2005 )). A minimal model, used to address particle Lagrangian dynamics 



in highly turbulent suspensions, takes into account only a few ingredients: the Stoke s 
drag, the Auton added mass and the fluid acceleration term ( Babiano et all ( 2000( ): 
Calzavarini et al. ( 20086f )). This leads to: 

dt Pf + 2 pp \ Dt J 

where pf and pp are respectively the fluid and the particle density, v the fluid kinematic 
viscosity and a the radius of the particle, which is considered spherical. However, the size 
of the particle in the above equation is essentially virtual. Equation (jl.ip contains only a 
time-scale, namely the particle relaxation time Tp, which embodies a particle length-scale 
merely in combination with the kinematic viscosity of the flow and with the densities 
coefficients, i.e., Tp = a^{pf + 2 pp)/{9vpf). In practice, the drag term in equation (|l.ip 
performs a purely temporal filtering of the flow velocity fluctuations. 

It is the role of Faxen terms to account for the non-uniformity of the flow at the 
particle-size. Faxen forces represent necessary physical corrections when analyzing the 
behavior oi D > rj particles in turbulence. The Faxen theorem for the drag force on a 
moving sphere states the relation 



tjj^&-Ki>Pfa (y-^^ J u{x) dS - = 6Tnypfa{{u)s^ - v) 



(1.2) 



where the integral is over surface of the sphere and u(x) the nonhomog eneous steady 
motion of the fluid in the absence of the sphere. As later shown bv Gatignoll ( 19831 ). 



Faxen force corrections via sphere volume averages should also be included on the inertial 
hydrodynamic forces acting on the sphere. In particular the expression for the fluid 
acceleration and added mass force becomes: 

where similarly as above (. . denotes the volume average over the spherical particle. 
Putting together the two force contributions of Eqs. (|1.2p - (ll.3p into an equation of motion 
for a sphere, (4/3)7ra'^Pp dv/dt — + f^, and keeping into account the Auton added 
mass correction for finite Rep, i.e., du/dt Du/Dt, we obtain the phenomcnological 
Faxen-corrected equation of motion: 

-77 = [{-^)v. + — i{n)s^-v)). (1.4) 

dt Pf + 2 Pp \ dt J 

In the small particle limit, when u ~ i;, corrections can be approximated by Taylor 

expansion (u(x,t))5„ ^ n + ^V^u + 0(0^); - ^ (u -f f^V^u + ©(a^)) , 

therefore the first order Faxen correction accounts for the curvature of the unperturbed 
flow at the particle location. In a turbulent flow the correction term becomes important 
when a > r], with a weak Taylor-based Reynolds number Rex dependence. With the as- 
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sumption V^u ~ 15urms/A^, where A is the Taylor microscale, one can roughly estimate 
a^V^u ~ 15u™s {a/\f ^ u™^ (a/rjf^/lE/ Re\. 



2. Numerical implementation of particle model and tubulence DNS 

We adopt here a further approximation which allows efficient numerical computations 
of Eq. (|1.4p . Volume averages at particles' positions are substituted by local interpolations 
after filtering by a Gaussian envelope with standard deviation, tr, proportional to the 
particle radius. Gaussian convolutions are then efficiently computed in spectral space, 
and the Gaussian volume averaged field reads: 

K)G,y„(x) - VTT-}^,^ [G<,(k) u,(k)] , (2.1) 

where VTT^^^^ denotes a discrete inverse Fourier transform on a grid A^"^, Gcr(k) = 

exp(— (7^k^/2) is the Fourier transform of a unit volume Gaussian function of variance 
a and ?Ii(k) is the Fourier transform of a vector field (the material derivative of fiuid 
velocity in Eq. (|1.4p '). The surface average is obtained using the exact relation: 



which leads to: 



{u,)g,s. (x) - -DTT-^,^ \S„{\^) ?I,(k)J , (2.3) 

where: ^^(k) = (l — -Icr^k^) e^^"'^^^ . It can be shown that with the choice a = a/-\/5, 
the Gaussian convolution gives the right prefactors for the Faxen correction in the limit 
a ^ 0. Our simplified approach for the integration of p.4p (EC model) allows to track 
inertial particles in turbulent flows with minimal additional computational costs as com- 
pared to Eq. p.ip (PP model) : the fiuid acceleration and velocity fields are filtered once 
for every particle radius size, then the averaged flow at the particle positions are obtained 
through a tri-linear interpolation. We track particles via Eq. (|1.4[) in a stationary homo- 
geneous isotropic flow, generated by large-scale volume forcing on a cubic domain. The 
Navier-Stokes equation is discretized on a regular grid, integrated using a pseudo-spectral 
algorithm, and advanced in time with a 2"'^ order Adams-Bashford integrator. 

We have explored in a systematic way the two-dimensional parameter space [pp/p/, D /rj] 
in the range pp/pf G [0.1, 10] and D/r] € [2, 50] for a turbulent fiow at Rex = 180 (512^ 
collocation points). We tracked '-^ 2 ■ 10^ particles for a total of ~ 4 large-eddy turnover 
time in statistically stationary conditions. Lower resolution DNS at Re\ = 75 (128^) 
have been used to explore a larger parameter space and to study the differences between 
the point-particle model (jl.ip and the Faxen corrected model (|1.4p in the asymptote 
D ^ L (with L the turbulence integral-scale). The validation of the numerics have been 
performed through careful comparison with an independent code implementing the same 
algorithm for particles, but with different forcing scheme, temporal integration method 
(Verlct algorithm), and local interpolation scheme (tri-cubic algorithm). 



3. Phenomenology of Point-Particle and Faxen-Corrected models 

We compare the statistics of acceleration of particles tracked via the PP and EC 
equations. In the small particle limit {D/rj 0) the two model equations behave the 
same way and the particle trajectory becomes the one of a fluid tracer. The ensemble- 
average acceleration variance reaches the value (a^) {a'i) with the subscript / labeling 
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Figure 1. (left) The acceleration variance (a^) normalized by /3^{aj) versus the particle diameter 
as derived from the Faxen-Corrected model (solid lines/symbols), and from Point-Particle the 
model system (dashed lines/empty symbols). Density ratios shown are Pp/pf = 0.1, 1, 10, i.e., 
heavy (□), neutral (o) and light (A) particles, (right) Same as above for the acceleration flatness 
F{a) = {a'^) / {a?)"^ . Horizontal lines shows the flatness of the fluid acceleration F{af) and the 
flatness value for Gaussian distribution F{a) — 3. Data from simulations at Rex = 75. 

the fluid tracer acceleration. As the particle diameter is increased we notice important 
differences between the two models. In the PP model the drag term becomes negligible 
and one gets (a^) ~ /3^(ay), with /3 = 3pf/{pf + 2pp). In the FC model the volume 
average of the fluid acceleration Du/Dt reduces progressively the particle acceleration. 
This is illustrated in Figure [1] (left), where the particle acceleration variance (normalized 
by /3^(aj)) is shown for three cases: neutral buoyant, heavy {pp/pf = 10) and light 
{pp/pf = 0.1) particles. We note that the behavior of (a^) for particles whose diameter 
is roughly larger than lOry seems to be identical apart from the scaling factor /3^. 

Differences are also present in higher order moments: for this we focus on the flatness 
F{a) = {a^)/{a'^y. In the large D limit PP model gives the rather unphysical behavior 
F{a) ~ F{af ), that is to say large particles, irrespectively of their density, show the same 
level of intermittency as a fluid tracer. On the other hand the Faxen corrected equation 
gives asymptotically F(a) ~ 3, i.e., the Gaussian flatness value, meaning that acceleration 
of large particles independently of their mass density value has lost its intermittent 
character, see Fig. [1] (right). Furthermore, it is noticeable that above a certain critical 
value of the diameter the flatness of heavy/neutral and light particles reaches the same 
level: this suggests that also the PDFs may have very similar shapes. 



4. Comparison with experiments 

We study now how the FC model compares with the experimental observations listed 
in the introduction - recalling than none is captured by the PP model. 

4.1. Acceleration variance 
In flgure [2] the behavior of the one-component acceleration variance, normalized by 
the Heisenberg-Yaglom scaling, ao = (af )e~^/^i^^/^, is displayed. Although this way of 
normaliz ing the a c celera tion has a weak Reynolds number dependence (see Voth et~al\ 
12002); Bee et al\ (j2006[ )) we notice a very similar behav i or as compared to the exper- 
imental measur ements at Re\ = 160 by lOureshi et all (|2007t ) and with Rex = 970 
experiments by IVoth et al. I (|2002l) . In the inset of Figure [5] the same quantity but with 
a different normalization is shown. The particle acceleration variance there is divided 
by the second rn oment of fluid tracer acceleration (a^). The experimental data from 
Voth et all (|2002h can also be rescaled in the same way by dividing oq by the value for 
the smallest considered particle (which has size D ~ 1.44ry and essentially behaves as 




Figure 2. One component acceleration variance versus particle size. Acceleration is normal- 
ized by the Heisenberg-Yaglom relation, while the particle size is is normalized by the dissipa- 
tive scale. DNS resu lts have uncertainty of the orde r of the symbol si ze. Data EXP are from 
iQureshi etal\ (|2007l ). with EXP* measurement from IVoth et al\ (|2002l ') (Figure 32) - particles 
with density contrast Pp/pf ~ 1.06. Inset: (a^)/(a^) vs. D/rj from the same DNS and EXP* 
measurements. 



a fluid tracer). This alternate way of looking at the data renormalizes the weak Re\ 
dependence, providing a good agreement between the DNS and experiments even when 
comparing results with one order of magnitude difference in Re\. 

In a DNS one can estimate the relative weight of the terms contributing to the total 
acceleration: the drag and fluid acceleration terms, respectively: — {{u)s — v)/Tp and 
a"^ = (3{^)v- It is important to note that in the case of neutrally buoyant particles, 
one finds {a)rms — (a"*)rms with percent accuracy. It indicates that the observed effect - 
decrease of particle acceleration variance for increasing particle diameter - comes uniquely 
from volume averaging of fluid acceleration at the particle position. The drag contribution 
is sub-leading at all D values (from few percent up to 15% of total acceleration variance); 
it just contributes to compensate the correlations. Stated differently, one can say 

that the acceleration of a finite-size neutrally buoyant particle is essentially given by 
{Du./ Dt)v^ — (V • T -I- fe)v„ — ^(r • n)5^, where t is the stress tensor, n a unit 
norm vector pointing outward the sphere and fg the external large-scale forcing whose 
contribution (fe)v'„ ~ is negligible at the particle scale. One expects the situation to be 
different for particles whose density does not match that of the fiuid. 

Our simulations arc consistent with the ao ~ D^^l^ scaling which has been proposed on 
the basis of dimensional arguments rooted on Kolmogo rov 1941 tu r bulen ce phenom enol- 
ogy w ithout special assumptions on particle dynamics ( Voth tt al. ( 2002 ): [Qureshi et al 



( 2007[ )): however at Re\ — 180 the scale-separation is still too limited to observe a true 



scaling range. 

4.2. Acceleration probability density function 

The second quantity under study is the acceleration probability density function. Here, 
to cope with Re\ effects, one compares only the t wo most simila r data sets: the DNS 
at Rex — 180 and the experiment at Re\ = 160 lOureshi et al. I (|2007l) . Experiments 



have revealed a universal behavior for acceleration PDF normalized by (a^)^/^ in the 
size-range I? = 12 — 25?;. DNS instead shows a systematic difference in the their trend: 




Figure 3. (top) Com parison of probability density functions of acceleration normalized by its 
rms value, from EXP iQureshi et'al] (|2007l ') at Rex = 160 and DNS at Rex = 180. (bottom) 
One component acceleration flatness F{a) = {at) / {O'l.)^ versus the normalized particle diameter, 
D/r), from the same experiment and DNS at two different Reynolds numbers. 



larger particles have less intermittent acceleration statistics, see Figure [3] (top). However, 
the shape of the PDF in the limit of large particles D ~ SOry shows a good similarity. 
To better visualize differences, in Figure [3] (bottom), we show the flatness F{a) vs. 
particle diameter for DNS and experiments. As already observed, the FC model leads 
to decreasing intermittency for bigger neutral particles, and in the asymptotic limit 
{D — > L) to Gaussian distribution; also acceleration flatness is an increasing function 
of Re\. Qureshi and coworkers' experiment on the other hand shows a D-independent 
behavior around F{a) = 8.5. A further possible source of differences can be connected to 
the variations in the large scale properties of turbulent flows: Experimental tracks come 
from a decaying grid-generated turbulence, simulations instead uses volume large-scale 
forced flow in a cubic domain without mean flow. 



4.3. Acceleration time- correlation 

Finally, we consider the dynamics of the neutral particles. We study the normali z ed one - 
component correlation function, Coa(T) = {ai{t)ai{t + T))/{af). In IVolk et al. ( 2008a ) 
it has been noted that PP model can not account for the increasing autocorrelation 
for larger particles. This is understood from equation (jl.ll) : In the large D limit the 
drag term is negligible and the acceleration of a neutrally buoyant particle is dominated 
by the inertial term Du/Dt. Therefore the time-correlation of acceleration, Caa{T), is 
related to the temporal correlation of Du/Dt along the particle trajectory. Because in 
the large D limit v ^ u ( Babiano et al. ( 2000l )). one expects an acceleration correlation 
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Figure 4. Autocorrelation function of acceleration Caaij) for neutral particles (pp = p/), 
with different sizes D = 2, 8, 16, 32r; and a tracer particle. Inset: integral acceleration time: 
Ti = Caa{r) dr with To the zero-crossing time Caa{To) = 0, versus particle diameter. 
Symbols: (□, o) for DNS at Rex = (75, 180), data (+) from experiments at Rex ~ 160. 



time which is equal or even shorter than the one of a fluid tracer. This is conflrmed by our 
numerics based on the PP equation (|l.ip . In the PC model instead, the averaged quantity 
{Du/Dt)v^ dominates the particle's acceleration and so its time correlation Gaai^)- In 
Figure |4] we show that simulations based on Eq. (11.41) display increasing correlation 
time for bigger particles, as observed in experiment ( Volk etal\ (2008a), although at 
much larger Re\ values). A detailed comparison of Caa{T) curves coming from DNSs 
with experiments by Qureshi and coworkers is at present not possible, because of limited 
statistics. Therefore, we examine integral quantities such as an integral acceleration time, 
T/. Since by kinematic constraint the time integral of Caaij) for a small tracer is zero, we 
define T/ as the integral over time of the positive part of Caa{T)\ this choice proves to be 
stable in the experiments and weakly dependent on the unavoidable (gaussian) smoothing 
of noisy data sets (see fVolk et al\ (|20086l )). The result of this analysis is reported in figure 
21 (inset). The order of magnitude ofTj/Trj, which is very near unity, as well its increasing 
trend with D qualitatively confirms the prediction of the PC model at similar Reynolds 
number. Using DNS results, it is also interesting to note that this time decreases with 
increasing Reynolds number. 



5. Discussion of results and conclusions 

We have investigated the origin of several experimental observations concerning neutrally- 
buoyant finite-size particle acceleration in turbulent flows and shown the relevance of 
Paxen corrections. Paxen terms account for inhomogeneities in the fluid flow at the spatial 
extension of the particle. They act as spatial coarse-graining of the surrounding turbulent 
flow, in contrast with the drag term which performs a temporal filtering. Numerically, the 
spatial average is efficiently implemented via Gaussian filtering in spectral space. Com- 
paring with experimental measurements, the main achievements of the Paxen-corrected 
model are: (i) prediction of the reduction of acceleration fluctuations at increasing the 
particle size; (ii) prediction of the increasing of acceleration time correlation at increasing 
the particle size. Both effects originate from the volume average of the fluid acceleration 
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term, or in other word from the surface average of the stress tensor of the unperturbed 
flow. While giving the correct trend, the FC model does not solve the puzzling point of 
invariant PDF with particle size, observed bv lQureshi et al\ (|2007l ). 

The Faxen corrected model marks a substantial improvement in the statistical descrip- 
tion of realistic turbulent particle suspensions. We emphasize that none of the observed 
trends in the acceleration of neutrally buoyant particles can be captured by purely lo- 
cal models, as e.g. the point-particle one in eq. (jl.ip . Faxen corrections are of highest 
relevance in the case of neutrally buoyant particles, because it is the case for which the 
slip velocity {vs = \v — (u)^!) is the smallest (as compared to pp ^ pf particles) and 
therefore where drag, history , lift have the least impact on the net force. In our case we 
observe that when increasing the size of particles, the PDFs of slip velocity normalized 
by the fluid velocity rms value {vg/urms) change from a sharp delta-like shape (for trac- 
ers) to larger distributions approaching a Gaussian (for large particles) . A size-dependent 
slip velocity for neutrally buoyant tracers in chaotic flows has been reported recently in 



Ouellette et al\ (j2008( ) : Faxen corrections to the added mass should be significant in that 



case too. We also observe that the particle-Reynolds number Rcp measured in our simu- 
lations attains values in the order of the hundreds, hen ce a more accurate description of 
the drag coefficient (see for instance Clift et all ( 1978h ) may be important particularly 
for a faithful reproduction of the far tails of the acceleration PDF. 
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